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INTRODUCTION 


The main tasks outlined in our original proposal to NASA "A Climatically- 
Derived Global Soil Moisture Data Set for Use in the GLAS Atmospheric Circulation 
Model Seasonal Cycle Experiment" included: 

1. The estimation of monthly, terrestrial soil moisture (w) fields 
at a resolution of 1° of latitude by 1° of longitude. 

2. The estimation and mapping (for verification purposes) of monthly 
potential evapotranspiration (E° m ) , actual evapotranspiration (Eg!), 
water surplus (Aw), air temperature (T m ) and precipitation (P m ). 

These fields also were to have a 1° x 1° resolution. 

3. The evaluation of Thornthwaite ' s E°m function in areas where it is 
thought to yield inaccurate estimates, e.g., in monsoonal and arctic 
environments. At the same time, the function which represents the 
resistance of the siil-plant system to E° m (8) was to be investigated. 

4. The final phase of the research was to be a comparison between the 
GLAS CCM-generated water balance and our Thornthwaite-based water 
balance . 

During the course of our research, the focus of our efforts evolved in response 
to unanticipated findings as well as to discussions with Yale Mintz and J. Shukla 
(the grant monitor) concerning NASA's changing requirements for soil moisture 
data sets. These considerations ultimately caused us to abandon any effort to 
complete Tasi. 4, for example, while we undertook a number of additional efforts. 
These additional tasks included the development of a global snow water equiva- 
lent data set as well as an algorithm which could efficiently and accurately 
interpolate from irregularly spaced station data to the nodes of a regular 
GCM-compatible lattice. Throughout our work, however, our overriding goal of 
estimating realistic monthly soil moisture values for the globe remained un- 
altered. Below we briefly describe our grant-related accomplisnments . Much 
of this work also was discussed by Willmott (1982) at NASA's Climate Science 
Review as well as in previous progress reports. 
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DATA SET COMPILATION 

Using average monthly air temperature and precipitation data previously 
compiled by C. W. Thornthvaite and associates, we created a digitally encoded 
global climate data base. The Thornthvaite data were edited and augmented 
such that our final data base contains 13,332 station records of monthly air 
temperature and precipitation as well as ancillary station information such as 
elevation, longitude and latitude. Willmott et al . (1981) published the data 
and a magnetic tape copy of the data plus documentation were made available to 
NASA's Goddard Laboratory for Atmospheric Science. The details concerning 
these data were not only described by Willmott et al . (1981) but they were dis- 
cussed in the second and third progress reports. 


ALGORITHM DEVELOPMENT 

When our proposal was written, neither ourselves or NASA personnel were 
aware that snow could not be adequately treated by our proposed water budget 
procedure. We further did not anticipate that adequate means were not available 

to interpolate values from the irregularly spaced st? cions to the grid points 

of a GCM-compatible lattice. Mapping procedures that were available, including 
the one we had proposed to employ (SYMAP), produced an unacceptable interpolation 
error because the interpolation computations were performed in a projected, 
Cartesian space. Our solutions to these problems, i.e., the development of 
new algorithms to include a snow budget in our water budget calculations and 
to interpolate and contour data in spherical space, are discussed by Willmott 

et at ■ (1984a) and Willmott et al . (1984b) (Appendix 1 and 2) as well as in 


the first, second and third progress reports. 
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EVALUATION OF THE THORNTHWAITE 
WATER BUDGET 


Evaluations of the Thornthwaite water budget were made by comparing (1) 
iysimeter-derived E 0 ^ estimates with E° m values predicted by the Thornthwaite 
model and (2) measured runoff, runoff predicted by one of Lattau's formulations 
and runoff derived from the Thornthwaite water balance (Appendix 3). The 
withdrawal function (@) was also investigated through sensitivity analysis 
(Willmott and Rowe, 1981) and evaluation with a modicum of published measure- 
ments . 

Comparisons between observed (Iysimeter-derived) E° m and Thornthwaite- 
prodicted E° m were inconclusive although scientific interpretation of results 
suggests that Thornthwaite 's E° ffi values are more representative of regional 
E° m than comoonly believed. Aspects of these comparisons are given by Willmott 
(1984) and they also are discussed in the third progress report. Evaluations 
of runoff predictions using measurements and Lettau's model were made by the 
Co-investigator, Richard Field, and his findings were similarly inconclusive 
(see Appendix 3). 

Plant-soil conductance to evapotranspiration is normally specified as a 
function of the ratio of soil moisture to field capacity (w/w^). This function 
is often referred to as the 'beta function' where f$ * f (w/w^). Willmott and 
Rowe (1981) found that the African seasonal soil moisture cycle, particularly 
in sub-Saharan Africa, was highly sensitive to the selection of the 8 “faction 
and their results are briefly outlined in the first progress report. Further 
research indicated that a 8-function suggested by Nappo (1975), when used in 
conjunction with a field capacity of 150 mm, would give reasonable soil moisture 
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and evapotranspiration estimates for much of the globe. Nappo's (197b) 6- 
function was subsequently employed in the global water budget analysis presented 
by Willmott et al . (1984b). 

ESTIMATED MONTHLY CLIMATE FIELDS 

Using the above-described air temperature and precipitation data set, and 
our water balance and mapping procedures, we derived monthly GCM-compatible 
fields of air temperature (T m ), precipitation (P m ), potential evapotranspiration 
(E° m ), actual evapotranspiration (Ejj), snow cover water equivalent (w 3 m ) and 
soil moisture (w^) at each of the 13,332 station locations. Values of these 
variables then were interpolated to the nodes of a 1° of latitude by 1° of 

g 

longitude lattice and visual verification (maps) of T m , P m , w m and w m were 
presented in the fourth progress report. A magnetic tape copy of the monthly 
fields of T m , P m , E 0 0 , Eu,, w m and w m was made available to NASA's Laboratory 
of Atmospheric Science. These data sets were further evaluated and described 
by Willmott et al . (1984b) and their creation marks the successful completion 
of the main task associated with our grant. An additional global water 
balance was computed from 899 of the World. Weather Records stations (at J. 
Shukla's request) and those fields also were magnetically encoded and given to 
the Laboratory of Atmospheric Science. This six-yesr (1974-1979) calculation 
is described in the third progress report. 


ONGOING RESEARCH 

Our research into the global water budget in part suggested that increased 
levels of biophysical realism should be incorporated into GCM-representations 
of the near-surface hydrology. The Principal Investigator with Yale Mintz and 
Piers Sellers of NASA subsequently requested and received funding from NSF and 



NASA to develop a model biosphere for use in GCM's (Mints etal. 1983). We 


are about six months into this research and we will likely use the above-* 
described data sets in the evaluation of our proposed biosphere model; that 
is, to see if our biosphere model improves GCM-predictions of near-surface 
hydrological cycle. 
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Algorithms for point interpolation and contouring on the surface of the sphere 
and in Cartesian two-space are developed from Shepard’s (1968) well-known, local- 
search method. These mapping procedures then are used to investigate the errors 
which appear on small-scale climate maps as a result of the all-too-common practice 
of interpolating — from irregularly spaced data points to the nodes of a regular lattice 
— and contouring in Cartesian two-space. Using mean annual air temperatures 
drawn from 100 irregularly spaced weather stations, the annual air temperature field 
over the western half of the northern hemisphere is estimated both on the sphere — 
assumed to be correct — and in Cartesian two-space. When the spherically- and 
Cartesian-approximated air temperature fields are mapped and compared, the magni- 
tudes (as large as 5 * to 10 ' C) and distribution of the errors associated with the latter 
approach become apparent. 

Key Words 


Climate, mapping, interpolation, contouring 
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Within climatology, automated point interpolation and isoline construction have 
widely replaced the relatively subjective methods of their draftsman predecessors. 
Cartographers have provided us with numerous digital algorithms for solving the typi- 
cal problems associated with interpolating from irregularly spaced data points in 
Cartesian two-space to the nodes of a regular lattice and isoline placement among the 
estimated grid point vaues (Peucker 1980; Rhind 1975; Morrison 1974; Marble 1981). 
It is possible, but uncommon in climatology, to draw isolines directly from the irregu- 
larly spaced data, ostensibly because the grid point values frequently are required for 
subsequent analyses. 

Another source of potentially useful interpolation and contouring procedures has 
been the atmospheric sciences where a wide variety of methods have been developed 
under the rubric of “objective analysis” (Panofsky 19^9; Cressman 1959; Barnes 1964; 
Fritsch 1971; Gandin 1963; Schlatter, et al. 1976; Wahba and Wendelberger 1980). In 
many ways, the development of objective analysis parallels the advancements in car- 
tographic interpolation and contouring except that the impetus for and direction of 
objective analytic research largely have been determined by the requirements of the 
weather forecasting community. Objective analysis, as a result, characteristically 
relies upon the typically smooth nature of upper air data as well as upon multivariate 
and/or autocorrelati^e relationships beyond those contained in the irregularly spaced 
scalar field associated with a single variable. Many objective functions, for these rea- 
sons, may be inappropriate for mapping the characteristically uneven scalar fields 
associated with near-surface climate although the objective analysis literature con- 
tains a number of innovative approaches to interpolation (e.g., Wahba 1979, 1981; 
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Wahba and Wendelberger 1980). 

It is important to ;.ote that most objective functions, similar to interpolation and 
contouring procedures developed by cartographers, tacitly assume that the data 
values have been projected into Cartesian two-space (e.g., Wahba and Wendelberger 
1980). Only in a minority of instances have the true relationships between the data 
and grid points on the surface of the Earth or sphere been preserved and explicitly 
incorporated into an objective algorithm (e.g., Wahba 1979, 1981). 

Point interpolation and isoline construction in climatology characteristically 
begin with the collection of data values associated with an irregularly distributed set 
of points on the surface of the Earth. Far most climatological purposes, the Earth 
can be assumed to Li a perfect sphere, with a trivial loss of accuracy, and each data 
observation (subscript i) of the variable of interest, z,(X,0), can be uniquely located 
by simple longitude (X) and latitude (0) coordinates. Once the data points have been 
compiled, they usually are projected into a Cartesian two-space so that each r t (X,0) 
becomes r,(x,y), and one of the readily available Cartesian-based interpolation and 
contouring algorithms (e.g., Marble 1981) may be applied to the data. Many times 
the projection is indirectly accomplished by the estimation of x and y from an exist- 
ing map or by the 3cded, but otherwise unaltered, assignment of X and <p r and y, 
respectively. Estimates of Zj{x,y) then are made at the nodes (subscript j) of a regu- 
lar lattice by point interpolation, which is ordinarily followed by the lacing of isolines 
through the lattice of interpolated values. During contour lacing, the points — grid 
points in this case — once again are assumed to be correctly related by the projected 
Cartesian geometry. Assumed Cartes, an relationships between the projected points, 
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in. other words, underlies th 5 interpolation of grid point values and often contributes 
twice to the estimation of isoline locations. 

When mapping climate fields that extend over large areas of the Earth’s surface, 
however, it is well known that spherical distance and/or directional relationships can- 
not be accurately preserved in the translation to Cartesian two-space. Subsequent 
interpolation and contouring — based on the projected relationships — will neces- 
sarily be in error. The magnitude of the errors, of course, will completely depend on 
the map projection selected, the distribution of the data and grid points, and the pro- 
perties cf the interpolation and contouring algorithm. It will be shown that this error 
can be marked. 

In order to avoid the propagation of such errors in small-scale mapping, climatol- 
ogists should perform their interpolation and contouring on the surface of the sphere, 
i.e., in “spherical space”. That is, the interpolation and contouring processes should 
depend only upon the spherical geometry that relates the grid and data points on the 
Earth’s surface. Only after the grid point values and contour positions have been 
determined on the sphere should they be projected onto a map. 

Even though most climatological papers do not describe or even cite the interpo- 
lation and/or contouring methods from which their climate maps were generated, 
errors that are derivative from the Cartesian-based interpolations are commonplace 
within the literature. Any of the numerous SYMAP (Shepard 1968J, IMSL (Akima 
1978) or SURFACE II (Sampson 1978) computed contour maps, for example, assume 
a Cartesian space. On any given small- scale climate map, particularly when the algo- 
rithm is not described, the nature and degree of error may be difficult to ascertain In 
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cases where the mapped area extends to one or more of the poles or over all 360 
degrees of longitude, however, two symptoms can easily be recognized. Cylindrical 
projections which have as their edges the north or south poles and the dateline, for 
instance, typically exhibit 1) multiple isolines crossing the poles and/or 2) contours 
which do not meet at the dateline (e.g., Schlatter, e£ al. 1876; Halem, et al. 1982; 
Rivin and Kulikov 1982). Suffice it to say that consideration of interpolation and/or 
contouring errors produced by the projection of data points into Cartesian space prior 
to grid point interpolation and contouring largely has been overlooked by climatolo- 
gists. 

In order to examine the degree to which Cartesian-based interpolation and con- 
touring can ill-effect the small-scale mapping and analysis of climaij fields, we present 
the results of two sets of map-analyses based upon two modified versions of Shepard’s 
(1968) interpolation algorithm. These procedures also were augmented by two con- 
tour lacing routines in order to permit the plotting of isolines. Shepard’s (1968) func- 
tion was selected as a benchmark procedure because of its wide-spread use (e.g., by 
1980 — in the form of SYMAP — it was implemented at well over 500 organizations) 
and because it represents a class of point interpolation algorithms that have an intui- 
tive appeal (Shepard 1968) and are “relatively robust’’ (Rhind 1975:299). 

Our first implementation of Shepard’s method is quite similar to that of SYMAP 
in as much as all interpolations and isoline determinations assume that the data and 
grid points are correctly related in Cartesian two-space. Our second, “spherically- 
based” version computes all interpolations and isoline positions from the actual data 
and grid point locations on the surface of the sphere. Longitude and latitude coordi- 
nates associated with each spherically-derived isoline are then projected along with all 
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other substantive elements of the map onto a Cartesian two-space. 

Following descriptions of these two algorithms — with an emphasis on the 
spherically-based procedures which have not been discussed previously — we describe 
results obtained from comparative mappings of an average annual air temperature 
field which was sampled from a recently compiled world climate data set (Willmott, et 
al. 1981). Our two implementations of Shepard’s (1968) function are described in the 
next section of this paper followed by an examination of two air temperature maps 
which were derived from combinations of the above-mentioned air temperature field, 
two well-known cylindrical projections and our two interpolation and contouring algo- 
rithms. 

Grid Point Interpolation 

Shepard’s (1968) interpolation function estimates a value at each node of a 
predetermined lattice from a small number of nearby data points. Developed for 
interpolation within Cartesian two-space, the procedure accounts for the distance and 
directional relationships between neighboring data and grid points (Figure 1). The 
algorithm also has a limited extrapolation capability that permits grid points to take 
on values outside the range of the data. Although there are a few important differ- 
ences between the algorithms described here and Shepard’s (1968) original function, 
the essence of his function has been maintained. 

Based upon a simple gravity hypothesis and other considerations (Shepard 1968, 
1983), the interpolation procedure requires that the value of each nearby data point 
influence the estimation of the associated grid point value by an amount proportional 
to their distance from one another. Weights are ascribed to three categories of 
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di stance according to 

r i /3<d i ' t <r j (1) 

d ,.t> r , 

where r ] is the search radius (defined below) and is the distance from grid point j 
to nearby data point k. The subscript k (and later /) is used to reference o^e of those 
data points near j which belongs to the set of data points that influence the interpo- 
lated value at j. Following Shepard (1968), the search radius is initially defined as a 
constant, i.e., as the radius of a circle (on either the plane or the sphere) which con- 
tains seven data points — on the average. In cases where less than four or more than 
10 data points fail within r ; of j, r ; is adjusted. When the former occurs, r } is set 
equal to the distance of the fifth closest data point while, in the latter case, r ; - is re- 
defined as the distance to the eleventh nearest data point. Four or 10 data points 
then will lie within the circle (area) of influence represented by the adjusted search 
radius. Each grid point value thus represents a constant area of influen t — except 
when fewer than four or more than 10 data points fall within the initial search radius 
— and from four to 10 data values. 

For each grid point, in other words, two sets of n ; (4<n ; <10) nearby data points 
are selected; one set on the basis of their Cartesian distances from j in the projected 
two- space and the other on the basis of their distances from j on the surface of the 
sphere. On the sphere (Figure la), d* k — the great circle distance between points j 



and k — is obtained from 
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cos d' k = sin 0 ; sin 0 t +cos ^ ; cos 0 t co: — \ k ) (2a) 

whereas Cartesian distance (Figure lb) is computed from the well-known rectangular 
distance formula 

d /,* = [(**— *;)*+(»*— » y )*]- 0 " 5 ^ 2b ) 

It should be noted that the superscripts c and s, respectively, are used to distinguish 
Cartesian- from spherically-derived distances and angles. When a distance or angle 
appears without a superscript, it represents either system of geometry. Each set of n ; 
nearest neighbors then is sorted in ascending order which provides that d ;1 and d J<Bj 
are associated with the nearby data points of minimum and maximum distance from 
j, respectively. 

Once ( S k , k—l t rij) have been computed (1), Shepard’s function corrects for the 
“directional isolation” of each data point relative to all the other nearby data points. 
Although a few authors (Morrison 1974) have suggested that such adjustments may 
be inconsequential or even deleterious to the accurate reproduction of a known 
functionally-derived surface, such a correction makes theoretical sense (Shepard 1983). 
Directional isolation of each nearby data point with regard to j is computed from 

[l -cos^fc.ol, (3) 

i=i 1 1 

where 9j{k,l) is the angular sepruration of nearby data points k and l when the vertex 
of the angle is defined as grid point j (Figure 1). The angular solitude of a data point 
with respect to j — on the sphere — is derived from 

cos *;(fc,0=(cos dh - cos d * tk cos <*/,,)( sin d^sin d' ti ) 


(4a) 
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For points projected into Cartesian space, the separation is obtained from the defini- 
tion of the angle between two vectors; that is, 


cos (*,/)= 


{x k -x j ){x l —Xj)+{y k -yj){y l —y j ) 




l^k 


(4b) 


Thus, data points that have a small angular separation contribute less individually 
than points which have a large angular separation. With the directional isolation cf 
data point k known, a combined set of weights is calculated from 


*t=S k 2 1 + 


TJ E 

i—i 


l^k 


(5) 


which limits the maximum influence of the directional isolation of A; to a doubling of 

the weight based on distance (Shepard 1983), i.e., when T k = S t . 

i — i 

In order to obtain non-sero gradients on the interpolated surface at data points, 
increments (Az*) are computed and added to the respective data point values (z*). 
The correction involves finding an average weighted gradient for each of the n ; data 
points within r ; of j, based upon the collective rates of change at the other data 
points within r } . On the surface of the sphere, the incremental corrections are 
obtained for each of the data points associated with j from 


A Zi 


-{ 


d(AzT 

dX 


d[(j,k)+ 


d(Az) 

90 


^}(M)|| U /( V 4* d ‘ tk ) } 


( 6 ) 


where 9{Az)/9X is the average partial derivative with respect to longitude, 9(Az)/90 
is the average partial derivative with respect to latitude, <*'(;,£)— (X ; — X t )cos <f> } 
and df{j,k)= <p ; — <f> k . The average partial derivatives are taken as 
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where dj[ (/,£)= (Xj — \ k )cos 4> k and dj(f»fr)= 0j — 0k, while a somewhat arbitrary 
adjustment is made in order to limit the influence of the increments to one-tenth the 
data range; that is, so |Az* j cannot exceed 0.1(max z, — min z f ) where again i refers 

to any data point on the map. Although it is not obvious, the limit of Jaz* j can be 

deduced from equation (6) by making use of the Cauchy-Schwarz inequality. The 
adjustment parameter (v) takes the form 


v = 0.1 (max z, — min z,) max 


21 — 0.5 


Forms of (6), (7) and (8) also are used for the Cartesian-based computations with 
dAz/dx, d x £ (j,k) ( dAz/dy, d‘(j,k), dfo, d^{l t k) t d${l,k) and d^ replacing dAz/dX, 
d{{j,k), dAz/d0, d'i(j,k), dj tk , d{{l,k), d'+{l,k) and respectively. 

With the evaluation of the above described functions, the value predicted at grid 
point j (2j) becomes 


E W k {z k +^z k )/ s ^ 
*=l k* =1 


d i.l> e 


m-‘ S 2 , 


d ,.i<‘- 


( 9 ) 
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where m is the number t * data points within e (defined below) of grid point j. That 
is to say, when one or more data points are sufficiently close to a grid point (i.e., 
within i of j ), the grid point takes on their average value; otherwise, part one of (9) 
defines the interpolation function. On the surface of the sphere, e is specified as the 
maximum of a function of 1) the grid cell longitudinal difference (AX) and the latitu- 
dinal range of the grid points , or 2) the grid cell latitudinal difference (A^); that is 


0.01max{0.5AX[co3 (max^)-f-cos (min^ ; )],A^} (maz^;)(min^-)^0.0 
0.01m&x{0.5AX[cos (maxl^y l)+l],A^} (max^jXnuii^)<0.0 


(10a) 


where xnax0 ; and min0 ; are the largest and smallest latitudes, respectively, found 
among all the grid points and maxl0 ; l is the grid latitude of greatest magnitude. In 
other words, when the grid points lie entirely within either the northern or southern 
hemisphere, e* is computed from part one of (10a) whereas part two of (10a) is 
evaluated if the grid spans the equator. For points in Cartesian two-space, 


c e =0.01max{ Ax , A j/ ) 


(10b) 


where Ax is the width and Ay is the height of a typical rectangular grid cell. Using 
these modified versions of Shepard’s (1968) algorithm, point interpolations from irreg- 
ularly spaced data points can be performed for any bounded rectangular or spherical 
lattice as well as for any spherical grid that spans the surface of the globe. 


Isoline Construction 

Once data values have been estimated at all nodes of a lattice, it becomes neces- 
sary to describe the gridded data field by fitting and subsequently plotting isolines. 
Cartesian-based methods of contour-lacing are well-known (e.g., see McCuilagh 1981) 
but a spherically-based lacing routine also was needed in order to position the isolines 
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on the surface of the sphere and, thereby, not contaminate the presentation of a 
spherically interpolated field with a projection-related lacing error. 

Our primary lacing unit is either a spherical (AX by A^) or Cartesian (Ax by 
Ay) grid cell which has its vertices defined by four grid points (Figure 2). Each grid 
cell is quartered by the definition of a grid cell “center” and the construction of 
“half-diagonals” which extend from each of the grid cell corners to the center. On 
the sphere, the center is located at the average latitude and longitude of the cell, 
while the half-diagonals are the arcs of great circles. In Cartesian two-space, however, 
the center and half-diagonals are defined by the intersection of the diagonals of a grid 
quadrilateral (cell) which has as its vertices four projected grid points. Within each 
spherical or Cartesian grid cell, four triangles are thus created. At the vertex which is 
shared by the four triangles (the cell center), a value is estimated as the arithmetic 
mean of the four corner values although this can be a source of error. When the cell- 
size is relatively small, however — as is the case with the maps discussed in the next 
section — the error is insignificant. On the other hand, when the lattice is course, the 
center point and its associated value should be more accurately determined. With the 
definition of a tnangular mesh and estimation of center print values, all desired iso- 
lines can be unambiguously positioned within the interpolated scalar field. 

From an initial point on or intersection with the edge of a first triangle, an isoline 
must pass through the triangle and intersect one of the two opposite sides at a point 
determined by the relative magnitudes of the values associated with the end-points of 
the intersected edge. The coordinates of each intersection are recorded. As the 
second intersection positions the isoline on an edge of a triangle adjacent to the first, 
the process may be repeated for the second triangle such that the coordinates of a 
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third intersection are stored. Then a third triangle adjacent to the second may be 
evaluated and so on until the isoline either crosses a map boundary or closes itself. 
By successively leproducing these calculations for the desired number of contour lev- 
els, a complete isoline representation of the gridded data field can be computed and 
stored 

From the series of intersection coordinates associated with the desired isolines, 
Cartesian-derived contour lines merely have to be scaled and plotted — along with 
such supportive information as a similarly projected graticule and/or land-area out- 
lines — in order to produce a map. The coordinate points associated with isolines 
that were constructed on the sphere, however, must be projected — in the same way 
as any background information — prior to scaling and plotting. 

A Comparison of Spherically- and Cartesian-Derived Annual Air Temperature 
Maps 

In order to illustrate the degree and nature of the error which can be produced on 
small-scale climate maps by projecting the data points prior to interpolation to a reg- 
ular lattice and contouring, annual average air temperature (*C) in the Northern 
hemisphere was sampled and alternately mapped by the above-described sphericaily- 
and Cartesian-based interpolation and contouring procedures. Using the data set 
compiled by Willmott et al. (1981), 100 stations were randomly drawn from the sub- 
set of stations bounded by — 170* <X< — 50* and 2*<^<90‘ which resulted in a 
markedly uneven distribution of station locations (Figures 3 and 4). Two map projec- 
tions — Lambert’s equal-area and Miller’s — representing the class of cylindrical map 
projections, were used to examine projection-related dissimilarities that can occur 
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bctween 1) isotherms that were estimated on the surface of the sphere prior to projec- 
tion and 2) isotherms that were determined in Cartesian two-space from a previously 
projected set of data points. Cylindrical projections were examined because they 
comprise the largest group of projections used to portray meso- and macroscale cli- 
mate Helds. Projection-related translations were performed by an experimental ver- 
sion of MAPRO (Kansas Geological Survey 1981) and the land-area outlines were 
taken from WORLDDATABANK I (Central Intelligence Agency 1972). 

Interpolation began with the superposition of a 4 * of latitude by 5 * of longitude 
lattice over the area to be mapped. A grid size of 4 * X 5 * was chosen because it is 
frequently used in global climate models. Grid point estimates of air temperature 
and, subsequently, isotherm positions first were computed for the spherical distance 
and directional relationships between the data and grid points. The X and <t> coordi- 
nates associated with each spherically-estimated isotherm and map were then two- 
times projected into Cartesian two-space, i.e., by Lambert’s equal-area and Miller's 
projections. Two comparable Caitesian-deri^ed isotherm maps next were constructed 
by producing two projected sets of data points — using Lambert’s and Miller’s pro- 
jections — from each of which Cartesian-based interpolations to a correspondingly 
projected (now rectangular) 4 * X 5 * grid then were made. The two Cartesian- 
estimated, grid point air temperature fields were each separately contoured within 
their respective projected Cartesian spaces. Once again, all computations were carried 
out by the modified versions of Shepard’s (1968) function and the isoline lacing algo- 
rithms described in the preceding sections. 

Before examining individual differences between spherically- and Cartesian- 
derived isotherm patterns, certain characteristics that the maps presented here have 
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in common — owing largely to assumptions embodied in Shepard’s (1968, 1983) algo- 
rithm — should be mentioned (Figures 3 and 4). Where the station network is sparse, 
both the spherical and Cartesian versions of Shepard’s function will extrapolate into 
the voids in such a way that the gradient of the temperature surface will diminish 
with distance from the nearest data points. This is why, for instance, the isotherms 
spread apart between Hawaii and Central America. In a few instances, where a large 
area contains but a single datum, the data point over-influences nearby and not-so- 
nearby grid points which in turn creates a flattening of the interpolated surface about 
the data point, e.g. around Hawaii. Near the edge of the flattening — as a conse- 
quence — the surface gradient must rapidly increase in order to account for the influ- 
ence of “new” data points that have come within the search radii of the local grid 
points. The steep north-south gradient between Hawaii and Alaska at approximately 
40 ' N latitude — evidenced by the convergence of the isotherms (Figures 3 and 4) — 
exemplifies this errant tendency. 

No attempt has been made to smooth the isotherms — in order to circumvent an 
extra computational expense — and, therefore, they have a slightly jagged appear- 
ance. At the scale these maps are presented, however, the lattice is relatively fine 
and, therefore, the maps’ appearance is not seriously compromised. The relatively 
fine texture of the grid further allows for a certain clarity of interpretation insofar as 
the err-n exhibited by the Cartesian-derived isolines can be ascribed to the interpola- 
tion process even though, generally speaking, isoline misplacements result from the 
compound of interpolation and contouring errors. 

When the relationships between the data and grid points on the sphere are 
preserved during the contouring process, the relative positions of isotherms remain 
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constant regardless of the map projection (Figures 3 and 4). In cases where the data 
points are projected prior to Cartesian-based interpolation and contouring, however, 
the locations of isotherms relative to each other and to the supplemental map infor- 
mation (e.g., land-area outlines and data station locations) will inconsistently change 
from one irregularly spaced station distribution to another as well as from one projec- 
tion to another. The inconsistency with which such errors are manifested from data 
set to data set and projection to projection is a result of synergism among the data 
point distribution, the pattern of the lattice and the inherent error associated with the 
map projection selected. That is, the error field on such a map may well be a unique 
result of the combined influences of the above three factors and, therefore, rather dif- 
ficult to evaluate. When isotherms also have been estimated on the sphere, however, 
they may be compared to the Cartesian-derived isotherms for any particular map of 
interest. 

Miller’s cylindrical projection — in combination with the data and grid point dis- 
tributions — produced the most striking differences between the spherically- and 
Cartesian-estimated isotherm locations (Figure 3). In the high latitudes — where the 
projection is highly exaggerated and there are few data points — there is a marked 
disparity between the two sets of isotherm patterns. When the corresponding 
spherically- and Cartesian-interpolated grid point values were subtracted from one 
another, the magnitudes of the differences in the higher latitudes (^>60*) were fre- 
quently in excess of 5 * or 10 * C. Tha Cartesian-derived isotherms generally fall north 
of their spherically-computed counterparts because the projection places the data 
points progressively further apart as the pole is approached. Another expression of 
the high-latitude error field is the Cartesian-derived —5' and —10 ‘C isotherms 
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which head northwest from about 70 * north latitude with the cooler of the two ulti- 
mately passing through the pole. Since it is well-known that annual isotherms ring 
the poles — as is evident among the spherically-derived isotherms — there is a very 
small possibility that the Cartesian-derived —5* and — *10‘C isotherms are correct. 
Such errors are common when Cartesian-based interpolations are made within a 
cylindrically-projected two-space although, in this instance, the errant tendency is 
somewhat intensified by Shepard’s constraint on the extrapolated surface gradient. 
Where the data points are most dense (e.g., in the United States) the agreement 
between the two isotherm systems is the best — although discrepancies still are 
apparent. 

On Lambert's cylindrical projection, where latitude is compressed as the pole is 
approached in order to preserve the equal area property, the errors are similar in mag- 
nitude to those associated with Miller’s projection (Figure 4). In this instance, how- 
ever, the data points are placed increasingly near one another as the pole is 
approached which, in turn, creates an unrealistically high local variance in air tem- 
perature zi the higher latitudes. As the Cartesian-based interpolation and contouring 
procedure treats the pole a. a lin'; on such cylindrical projections, the abnormally high 
variance in temperature causes large fluctuations in the interpolated field and its 
isothermal representation. Most notably, this can erroneously cause multiple isoth- 
erms to pass through the pole. Once again, within the area whe’-e the distribution of 
data points is the most dense, the Cartesian- and spherically-derived isotherms agree, 
but not perfectly. 

Symptoms such as these are typical of small-scale, Cartesian- based maps which 
employ cylindrical projections, especially wb*" such a map surface includes the higher 
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latitudes. As mentioned above, when a cylirdricsJly-projected map encompasses the 
full 360 ' range of longitude, it also is quite common for east- west trending isolines not 
to connect at the common longitude represented by the map’s vertical edges. It 
should be emphasized, however, that such errors also exist on c’l small-scale, 
Cartesian-derived isoline maps although their identification often is difficult because 
1) the map does not extend to at least one of the poles or over 360’ of longitude, ?,) 
the interpolation and/or contouring process may be constrained near the map borers 
in order to make the map appear correct, or 3) too little information about the inter- 
polation and contouring methods is provided for the map-reader to make a meaning- 
ful evaluation. For these reasons, we recommend that interpolation ’ o a regular lat- 
tice and contouring should take place on the surface of the sphere after which the iso- 
lines and other map information may be projected. 

Summary and Conclusion 

A sensitivity study of the errors which can occur on small-scale climate maps as a 
result of the common practice of projecting the data points prior to interpolation i nd 
contouring has beer presented. Working from Shepard’s (1968, 1983) well-known, 
local-search interpolation function, two algorithms that perform the interpolation and 
contouring process 1) on the surface of the sphere and 2) in Cartesian two-space are 
developed and described. In order to illustrate the natice and magnitude of the 
differences that can occur between spherically- and Cartesian-derived interpolations 
and contours, mean annual air temperature over the western half of the northern 
hemisphere — represented by 100 randomly selected stations — was alternately 
mapped by the two algorithms on *wo well-known map objections. 
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• When the spherically-derived isotherm maps — assumed to be correct — were 
superimposed over the Cartesian-interpolated and contoured maps, the error fields 
manifested by the latter approach became apparent. With regard to our air tempera- 
ture examples, local errors as large as 5 ’ to 10 * C commonly appeared in those regions 
that contained few data stations or had large projection distortion. When a 
Cartesian-derived, small-scale climate map encompasses most of +>, e globe, however, 
two characteristic errors often can be observed without the benefit of a spherically- 
derived standard. On such maps, it is common to find more than one iso.'ne passing 
through a pole or east-west trending isolines not to meet at the common longitude 
portrayed by the vertical edges of the map. 

Our findings strongly suggest that the interpolation to grid points and subse- 
quent contour lacing on small-scale climate maps should, in virtually all cases, be car- 
ried out on the surface of the sphere. To do otherwise is to assure that calculations 
performed on the interpolated grid point values or climatic inferences made from the 
isoline representations of the raw and gridded data fields will be in error. Since our 
spherically- and Cartesian-based version of Shepard’s algorithm and their accompany- 
ing contour-lacing procedures are logically identical — save that one set computes 
spherical distances and directions and the other Cartesian — our findings further sug- 
gest that comparable errors will result from the use of any Cartesian-based interpola- 


tion and contouring method. 
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Figure Captions 

Figure 1. Schematic representations of an irregularly distributed set of data points 
(C’*0 and of grid nodes (D’s) a) on the sphere and b) projected into 
Cartesian two-space. Also illustrated are the corresponding distance 
(d* tk and dj tk ) and directional (8*(k t l) and 8'(k,l)) relationships between 

the data and grid points in spherical (superscript s) and Cartesian 
(superscript e) space. It should be noted that 8j(k,l)^6‘(k,l) or 
d i,k^lk both. 

Figure 2. Schematic diagrams depicting the placement of a) spherical and b) 
Cartesian grid cell centers (A’s) and the subdivision of each grid cell into 
four triangular sub-cells. A hypothetical isoline (dashed) has been laced 
through the triangular mesh — intersecting the edges of triangular sub- 
cells at points (O’s) determined by either great circle or linear interpola- 
tion between the end-point values of the intersected edges — in order to 
illustrate the contour placement procedure. 

Figure 3. Isothermal representations of mean annual shelter-height air tempera- 
ture ( ' C) over the western half of the northern hemisphere using Miller’s 
cylindrical projection. From 109 irregularly spaced data points (A’s), 
temperatures were interpolated to a regular 4 * x 5 * lattice and subse- 
quently contoured. The solid isotherms represent the interpolation and 
contouring process as performed on the surface of the sphere prior to 
Miner’s projection whereas the dashed isotherms depict the interpolation 
and contouring process as carried out in Cartesian two-space following 
Miller’s projection of the data points into two-space. 
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Figure 4. Isothermal representations of mean annual shelter-height air tempera- 
ture (*C) over the western half of the northern hemisphere using 
Lambert’s cylindrical equal-area projection. From 100 irregularly spaced 
data points (A’s), temperatures were interpolated to a regular 4’ X 5’ 
lattice and subsequently contoured. The solid isotherms represent the 
interpolation and contouring process as performed on the surface of the 
sphere pri?r to Lambert’s projection whereas the dashed isotherms dep- 
ict the interpolation and contouring process as carried out in Cartesian 
two-space following Lambert’s projection of the data points into two- 
space. 
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ABSTRACT 

A zeroth-order Approximation of the large-scale spatial and seasonal vari- 
ability associated with terrestrial snow cover, soil moisture and evapotranspi- 
ration is presented and evaluated. These surface climate fields were obtained 
th "Mgh (1) an empirically-based water budget analysis of 13,332 station records 
which contained average monthly air temperature and precipitation observations 
19fil) and (2) the spatial interpolation of the estimates to 
regular, general circulation model-compatible lattices (1° of latitude by 1° of 
longitude and 4° by 5°). The grid-point seasonal cycles associated with the 4° 
x 5° mesh then were harmonically decomposed and the first two harmonics evaluated. 
Seasonal snow cover exhibited the most coherent space-time distributions whereas 

* v 

the soil moisture cycle was the most variable. Evapotranspiration displayed 
moderate levels of space-time variability. 



INTRODUCTION 


Sensitivity experiments with atmospheric general circulation models (GCMs) 
repeatedly have demonstrated the profound influence which the near-surface, 
terrestrial water balance exerts on climate (Mints, 1982). Available soil 
moisture and surface water (e.g., snow cover) are prircipal agents in the 
determination of regional albedo and, perhaps more importantly, they largely 
control the partitioning of the surface enthalpy flux between latent and sen- 
sible heat. The degree to which surface enthalpy flux is comprised of latent 
or sensible heat has climatic significance for the following reasons: sensible 

heat transfer from the terrestrial surface warms the local boundary layer whereas 
latent heat is usually realized in the free atmosphere during the condensation 
process— frequently at considerable distance from where it is added to the 
boundary layer. Ensuing dissimilarities in the vertical, spatial and temporal 
distributions of these two sources of atmospheric heat make the general circulation 
of the atmosphere (i.e., the thermal ly^Lnduced large-scale motion) sensitive to 
the near-surface convective fluxes. The moisture associated with latent heat 
flux from the land surfaces^evapotranspiration — is additionally important because 
it represents the major source of moisture for terrestrial precipitation (Mintz, 
1982). 

Since most GCMs use a one- or two "bucket" representation of the near- 
surface water balance (Carson, 1981) and it has been demonstrated that this 
over-simplification .3 prone to error (Mintz et al., 1983: Sellers, 1981), 
improvements in GCM-representations of the near-surface water balance are 
prerequisite to a fuller understanding of atmospheric circulation and climate. 

More realistic analogs of the role that the terrestrial biosphere plays in the 
water balance promise to remove much of the uncertainty and inaccuracy 
associated with the open-bucket approach (Dickinson, 1983: Mintz et al., 1983). 
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Anothcr problem has been and continues to be the paucity of reliable, 
observed global water balance data which are available for the operational and 
scientific evaluation of GCM predictions of the seasonal water cycle as well as 
for use as sub-annual initial or boundary conditions by GCHs. Jaeger's (1983) 
discussion of global precipitation data sets well illustrates this problem. 
Observed climate fields have additional utility in that inductive evaluations 
may reveal patterns and provide insights which may lead to improved GCM-repre- 
sentations of the climate system. Sub-annual precipitation and air temperature 
fields (seasonal and monthly) have been observed and compiled (e.g., Moller, 
1951: Jaeger, 1976: Steinhauser, 1979: Willmott et al ., 1981) but large-scale, 
subannual fields of evaootranspiration, soil moisture and snow cover are not 
well-represented in the literature. Baumgartner (1981), for instance— in his 
review of the large-scale water balance literature — noted that "World wide 
water balances for shorter periods than a year are not available at the moment." 
Efforts to upgrade both GCM-characterizations of the near-surface water budget 
(e.g., Hints et al. , 1983) and the store of empirically-derived sub-annual 
water budget observations (e.g., Willmott et al . . 1981) are ongoing by the 
authors although this paper is restricted to the presentation of heretofore 
unavailable monthly water budget data and inferences that can be drawn from 
those data. 

Mid-monthly values of terrestrial soil moisture, snowpack water equivalent 
and evapotranspiration were derived from a global air temperature and precipi- 
tation data set (Willmott et al . , 1981) by means of an empirically-based water 
budget algorithm. Since the water budget method employed represents a sub- 
stantial revision of a previously published procedure (Willmott, 1977), it will 
be described in the next section of the paper with an emphasis on pertinent 
upgrades. Our estimates of the seasonal soil moisture, snow cover and evapo- 
transpiration cycles as they vary across the world's land surfaces are presented 
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and evaluated in third section of this paper and concluding remarks follow. 


CLIMATIC WATER BUDGET 


A simple water balance algorithm, which could be used with monthly air 
temperature and precipitation observations taken from any terrestrial location, 
was developed from a similar, existing model (Willmott, 1977). The procedure 
assumes that monthly evapotranspiration (E m ) can be adequately evaluated by a 
calculation of the form 

Ej| * 8 ]|E 0 ib so month - * 

where represents the integrated conductance of the biosphere (near-surface 
plant-soil-water-atmosphere system) to evapotranspiration and E° m is reference 
crop or potential evapotranspiration. Computations are made daily and summed 
in order to obtain monthly totals. 

Unadjusted reference crop evapotranspiration is first estimated according 
to Thornthwaite' s (1948) method; that is. 
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where E°ni is set equal to zero when (T <0), is the average shelter height 
air temperature (°C) for month m, I is the station heat index and a is an empiri 
cally-derived constant. The latter two terms are obtained from 


12N _1 Z (T /5) 1 ' 514 

i ra 

m » 1 


a - 6.75 x 10’ 7 I 3 - 7.71 x 10 _5 I 2 
+ 1.79 x 10" 2 I + 0.49 


and 


-4- 


wtiere 17 is Che number of months or days from which a representative station 
heat index My be calculated. It (N) is usually a multiple of 12 or 365 for 

I 

monthly or daily budgets, respectively. After E° m has been computed, it is <- 
adjusted for variable day and month length according to 

E° - E°' 1(0 /30) (0 ,/12) J mm month' 1 

■ am c 

where is the length of month a (days) nd is the length of day (hours). 

Thornthwaite ' s (1948) estimate ot regional E° m occasionally has been 
criticized because it is not biopcysically-based (e.g., Lee, 1978) and it is 
thought to yield insufficiently accurate values of E 0 ,,. (Jensen, 1973). Such 
evaluations seem quite reasonable '.pon cursory examination but more careful 
investigation suggests that a our >e r ~ of the comsonly-leveled criticisms may be 
over-stated. It should also be kept in mind that Thornthwaite' s estimat ;s only 
require average air temperature data which are widely available and which represent 
fewer input requirements than most competitive methodologies (e.g., see Jensen, 
1973). 

By making E° a an empirical function of T^ , Thornthwaite (1948) was not 
ignoring the salient biophysics of which he was well-aware (Thornthwaite and 
Holzman, 1939) but rather he was attempting to make use of a weather variable 
that represented the integration of all the flux terms pertinent to the estima- 
tion of regional E° m . Too often, it is forgotten that E° m responds to warm and 
cool air advection as well as to irradiance and that the combined, temporally- 
integrated influences of these two fluxes on E° a covaries with the vicissitudes 
in T a . In other words, while the form of Thornthwaite r s model can only be 
partially justified on biophysical grounds and the magnitudes of his coefficients 
have virtually no biophysical interpretation, his selection of as the inde- 
pendent variable, a« opposed to net irradiance for instance, i±> wall-taken on 
both biophysical and practical grounds. In the end, however — Thornthwaite ' s 
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biophysical insights notwithstanding— it primarily is a mathematically-baaed 
model (Willmott, 1984a) which was pragmatically developed for predictive, not 
explanatory , purposes. As a consequence, the main criteria on which this formu- 
lation should be evaluated include: economy, ease of use and, most importantly, 
accuracy. 

Thornthwaite's model is often thought to produce inaccurate estimates of 
E° m —relative to many of the competing procedures — and this conclusion is 
based largely upon comparisons with lysimeter observations. Jensen (1973), for 
example, compared 18 frequently used methods of estimating E° m with monthly 
lysimeter observations drawn from ten sites and he concluded that Thornthwaite's 
model “did not produce estimates that agreed well with measured values," i.e. , 
the estimates were virtually all low. On the average, the root mean square 
error (BKSE) was 1.84 mm day - *. If the measured values are correct, this error 
would be unquestionably large; however, since lysimeters often over/estimate 
regional E° a , it is likely Chat Thornthwaite's values are closer to actual, 
regional E° B than commonly believed. In fact, if one assumes that the systematic 
difference between Thornthwaite ^-ia^and the lysimeter estimates of E° m can be 
attributed to the lysimeter-derived values, the RMSE expectation reduces to 
0.92 tm day - * (Willmott, 1984b) which is quite competative (Jensen, 1973). For 
the sake of brevity, the above discussion somewhat oversimplies an ongoing 
evaluation of Thornthwaite's method; nevertheless, these and other findings 
suggest that the Thornthwaite model provides more accurate estimates of regional 
E° m than is generally believed. 

In addition to the Thornthwaite estimate of E° m , the water budget calcu- 
lations require monthly precipitation information which is assumed to be in 
liquid form (P r m ) when (f > -1 ) ; otherwise, precipitation is treated as snowfall 

111 IQ 
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(P* a ). Estimates of E° a , P r m and P s a remain constant over the course of mon.:h 
m whereas the storage terms ( i.e., soil moisture^ stferwater equivalent! and 
actual evapo transpiration are computed daily in order to partially account for 
the nonlinear, time-dependent relationships between E° m , P 1 ^ , P* m and actual 
monthly evapotranspiration—Eg,. 

Terrestrial water storage is partitioned between a snow cover and a soil 
moisture store where the latter has a saturation or field capacity (w*) which 


must be specified prior to calculation. In this instance, w* was held constant 

I 

A 


at 150 an which is large for m soils that underlie grasses and other low 


vegetation, and small for oust forest-covered soils. Sensitivity comparisons 
between w* * 150 and more representative biome-specif ic values suggest^** 150 
yields an acceptable error response in nearly all biomes. Since E°,, , P r m and 
P®^ are taken to be constant over month m and the water budget is evaluated on 
approximately a daily basis (i.e., 30 times per month), the daily estimates of 
reference crop evapotranspi ration, rain- and snowfall are E 0 ^ ■ E° m /30 , P 1 ^ 
® P r n /30 and P*^ * P s a /30 where d refers to the day of hypothetical 30-day 


month m. 

Quasi-daily evaluation of the snow and soil moisture stores begins with 
the estimation of the water' contained in the snowpack. At the end of day d, 
the water content of the snowpack is 


w ®md * ^md-l + p8 md ' 

where w 3 ^.} is the water equivalent of the snowpack at the end of the previous 
(d-1) day and Mod is the snowmelt during day d. Daily snowmelt is estimated 
from 


mm day 


-1 


Mmd » 2.63 + 2.55 ♦ 0.0912 
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where "ad is constrained so that it canno'w be less than zero or greater than 
(t^md-l + P*^). Comparisons between the above equation and the data from 
which it was derived— representing three dissimilar drainage basins and 113 
observations (Anderson, 1973: Pysklywec et al., 1968: Storr, 1978) — suggest a 
moderately good fit, i.e., RMSE * 6.62 ssn day“*and r^ » 0.59. Evaporative 
demand from the surface or soil water surplus (Dj^j) subsequently is calculated 
as 

Dmd ■ IW + P r md “ E 0 ^ • mm day ~ l 


page is 
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When ®md is negative, it represents a demand and it constitutes a surplus when 
it is positive. The degree to which the biosphere facilitates evapotranspi- 
ration can be estimated according to 

1 - exp 0-6. 68 w^/w*) » D md < 0 


®md 


• D »d i 0 


where w^j.} is the soil moisture at the beginning of day d. Using data from 
Davies and Allen (1973), the constant (6.08) was translated from a best fit described 
by Nappo (1975). Available soil moisture at the conclusion of day d then becomes 


w md * w md-l + Sod D md * «■ 

If w^j exceeds w + , the difference is retained as a surplus (s^j) for day d 
and w^ is set equal to w^. The increase or decrease soil moisture during month 
m then can be defined as 

A«m " w m30 “ w m-l , 30. “» 


When the temperature and precipitation tioieseries (vectors) are continuous, 
span one or more complete years and each element of those vectors is a 
climatic normal or representative average; a water balance may be assumed. 


In such cases i an iterative solution to the above-described set of 
equations is possible where the iteration continues until the quasi-daily soil 
moisture and snow cover vectors (v and w s respectively) become virtually in- 
variant. Mote that each of these vectors has 301^ elements where Nm is the 
number of months over which the water balance is computed. Actual monthly evapo- 
transpiration then may be obtained from 


30 


30 


s 


E - P „ x I M - aw - l 
m m . _ , ad m , , md 

a * l d ■ 1 


mm month 
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and the monthly soil moisture deficit is 6^ * E° m - E m * Representative monthly 
soil moisture and snow cover values then are taken as those elements of w and 
w s that correspond to the fifteenth of each month evaluated. 


TERRESTRIAL SEASONAL WATER CYCLE 

Using the procedure described in the above section, annual water 

balances were computed for 13,332 stations from average monthly air 

temperature and precipitation data (Willmott et al ., 1981). At each station, the 

water balance was iteratively solved and the solution sets of monthly total 

evapotranspiration (E^ , mid-monthly soil moisture (w m ) and mid-monthly 

water equivalent of the snowpack (w s ) are evaluated. 

n 

Following their estimation, each irregularly-spaced climate field, e.g., 
the January evapotranspiration estimates (E}) associated with the 13,332 station 
locations, was used to spatially interpolate additional values at the 

O'? 

nodes of a 1° of latitude by l° A longitude lattice. The interpolations were <: 
performed by a method described by Willmott e£ al. (1984) with the exception 
that the number of "nearby" data points (stations) which influenced a grid point 
value was held constant at 10. This modification of Willmott al. (1984) 
has the effect of dampening the spatial 'aciance among grid points when 
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s tat ions are sparsely located and accentuating the local variability when 
stations are in close proximity to one another. In order to preserve map 
clarity, we report on a subset of the interpolated water budget fields that 
is associated with the nodes of a 4° of latitude by 5° of longitude grid — which 
is frequently used in GCM applications. 

Snow Water Equivalent 

Statistics that summarize the water balance-predicted seasonal progression 
of snow cover (in a water equivalent depth of nan) are generally consistent 
with climatological expectation and with available observations (e.g., Dewey 
and Heim, 1981), although no truly comparable data set exists to the authors' 
knowledge. The average quantity of water entrained in the snowpack ( V s ) 
increases with the magnitude of latitude, particularlyjin the Northern Hemisphere, 
although certain mesoscale departures are apparent (Figure la). Greater average 
snow depths, for example, occur within the continental climates of eastern North 

America owing to the synergy between relatively high precipitation rates in winter 

I 

and the A temperatures that persist in the absence of a marine influence. The accen 
tuating impact of orography and the lower air temperatures of the mid-latitude 
mountain ranges are evident, for instance, in the Rocky and Ural mountains. 

Owing to the general poleward decrease in land area and elevation in the Southern 
Hemisphere, this region — by contrasty has few locations with a seasonal snow 
cycle (Koliakov and Krenke, 1981) or an average snow water equivalent of more than 
10 mm (Figure la). Antarctica, of course, is the exception and the water balance 
estimates a variable but perennial snow-ice cover which, at virtually all locales, 
is well in excess of 500 tmn. 

Expressed as the standard deviation associated with the seasonal snow 
cover cycle (0 ), the spatial distribution of temporal variability also is 

consistent with climatological expectation (Figure lb). Where the snow cover 
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is seasonal (i.e., ic completely melts in the summer), a is colinear with 
9 s since v 8 a is bounded on the low end by zero. Moving northward, along the 
east coast of North America for example, the seasonal variance increases 

*3 

with latitude to approximately 5G°N where o ’begins to lessen because of the 
decline in winter precipitation and the reduction of the melt season length. 
Mountainous regions exhibit large seasonal variances as well as large average 
snow covers (Figure lb). Antarctica has relatively small variances that, 
once again, decrease with increases in the magnitude of latitude. 

Although the spatial distributions of the first and second moments asso- 
ciated with the seasonal snow cover cycle (Figure la and lb) weil-describe 
average snow cover and the degree to which snow cover varies seasonally, they 
do not provide information about the temporal progression of the seasonal 
cycle; for instance, the times of the year when extrema occur and the amplitudes 

' r 

associated with the main periodicities. These aspects of the seasonal cycle, 
however, can be described by means of the harmonic decomposition of the 
monthly timeseries. A location's seasonal snow cover cycle may be represented 

by 6 

S — s 3 

w ■ w + £ v, cos[ ir(km - ra' )/6] mm 

k « 1 K k 

where w 8 ^ is tne amplitude of the k-th harmonic, k is the frequency and m'^ is 
the phase angle in months. The parameters w 8 ^ and m'^ (k * 1, 6) then provide 
a complete description of the seasonal cycle. Since the higher-frequency 

. . . . . ~s 

harmonics typically explain local or trivial portions of o their interpre- 
tation is often of limited value. For this reason, we present only the 
parameters associated with the first two harmonics which, at the mesoscale, 
introduce an acceptable error-level. Over the continents of both the Northern 
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and Southern Hemispheres, the magnitude of the variance left unexplained by the 
first two harmonics [Cl - r^)o 8 ^]^*^ usually is not much greater than 10 mm 
(Figure lc). 

Phase angles associated with the first harmonic (k ■ 1) clearly represent 
the end of the accumulation season (winter) and the onset of melt (Figure Id). 

A clockwise shift in m'^ from January and February, near the snowline, to March 
and April in the higher latitudes is evident in both North Americaj^and Jiorthern ^ 
Asia. Regions such as eastern North America, that are little-influenced by 
maritime air in the winter, also exhibit an extended accumulation season. 

Local increases in both the water equivalent of the snowpack and the length of 
the accumulation season additionally appear in mountainous areas (e.g., the 
Sierra Nevada, the Pyrenees and the Himalayas). Consistent with the spatial 
distribution of a , the amplitudes associated with the first harmonic are 
greatest in the latitude band where winter, frontal snow most frequently fails. 
Similar patterns also appear in Antarctica but, of course, they are four or 
five months out of phase (Figure Id). 

Phase shifts (m'^/k) associated with the second snow cover harmonic 
(Figure le) are approximately equivalent to the phase angles of the first 
harmonic— near the equatorward limits of the seasonal snowpack. Angular de- 
partures of m' 2/2 from m' i , however, increase with the magnitude of latitude. 

The rotation of these departures is clockwise in the Northern Hemisphere and 
counterclockwise in the Southern Hemisphere. Second harmonic amplitudes also 
covary positively with the magnitude of latitude. In the lower latitudes^rhere 
m'i * m'2/2, the "direction" of m' 2 / 2 , relative to n'i , indicates that the 
temporal range of the snow cover is confined to a relatively short snow 
accumulation season. As tha magnitude of latitude increases, however, m'2/2 
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rotates in « clockwise precession from m r j in order to account for the lengthen- 
ing of thn snow accumulation season. A concomitant increase in the snow water 
equivalent (represented by the amplitude) confirms the increasing importance of 
this extension in the length of the snow cover season with latitude. Although 
most of the soatial variance in (m' 2/2 - m'j) can be explained by latitude 
alone, the melt process is regulated by the available energy and, therefore, 
the patterns are synoptically modified in some cases. Note the shortening of 
the snow cover season in northern Scandanavia relative to areas at the same 
latitude in eastern Siberia due to maritime warming. 

Soil Moisture 

Estimated mean soil moisture is much more spatially variable than average 
snow cover owing to its dependency on and interaction with a greater number of 
"independent" variables such as field capacity, biospheric conductance (B m ) 
and reference crop evapotranspiration (Figure 2a). Certain large scale patterns, 
however, are apparent and are consistent with well-accepted climatic principles. 
The humid eastern portion of North America, for example, has plentiful soil 
moisture throughout most of the year as manifested by averages well-above 
125 mm— assuming a field capacity of 150 mu. By contrast, the drier reaches of 
the western plains states as well as oi the Sierra Nevada average 8 25 mm. 

Orography along the northwest coast provides a thin ribbon of land with seasonally 
high precipitation (winter maximum) which is manifested in much higher soil 
moisture averages than the nearby mountain and valley regiona that lie to the 
east. 

South America too exhibits familiar patterns (Figure 2a). A very dry west 
coast (<l 25 mm) extends from northern Peru (near the equator) to central 

Chile and then southeast into Patagonia which lies in the rainshadow of the 
southern Andes mountains. The windward side of the southern Andes and the 
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coastal zone, on the other hand, are quite moist. In eastern Brazil, the 
tropical rainforest exhibits similarly high average soil moistures owing to 
excesses of precipitation but these extend over a much greater area. High 
average soil moistures along the northeast coastal reaches of South America 
as well as in southeastern Brazil also are correlated with high rainfall. 

Within Europe, average annual soil moisture is spatially variable although, 
at the regional scale, few areas can be classed as dry (Figure 2a). Most of 


f 

Europe exhibits plentiful precipitation which continually replenishes the soil 
moisture. Where eastern Europe grades into western and central Asia (i.e., 
beginning with the Caspian desert and the Ural mountains), precipitation is 
drastically reduced and the accompanying decreases in annual soil moisture are 
apparent over much of north-central Asia. Along the seasonally humid ear 
coast of Asia, however, generally high— but locally variable— levels of 
precipitation provide relatively high levels of soil moisture (>125 mm) in 
Japan, Korea and areas north. South and Southeast Asia receive even higher 
rates of precipitation and, therefore, their soils are usually near saturation. 
Much of this rainfall is convectional (e.g., in Sumatra, Java and Borneo) 
although in Monsoon Asia (e.g., in Burma, Bangladesh and Vietnam) plentiful 
precipitation is dynamically (i.e., mousoonally and orographically) induced 
during the suuner. 

Annual soil moisture patterns in Africa also spatially covary with 
precipitation (Figure 2a). Most striking, of course, is the Sahara Desert and 
circur -Sahara regions which together span northern Africa excepting the Atlas 
mountains and coastal areas of Morocco and northern Algeria. South of the 
Sahara, over much of central Africa as well as the Ivory Coast, the soil 
remains moist throughout the year ( >75 mm) principally owing to convective 


ORIGINAL PAGE IS 
OF POOR QUALITY 


-14- 

reinfall . Southwes tarn Africa, e.g., Che Kalahari desert, by contrast is c^uite 
dry throughout the year. 

Australia is very dry (S < 25) over most of its interior (Figure 2a). The 
northern, eastern and southern coastal regions represent virtually the only 
moist regions on the continent. 

Soil moisture variance associated with the seasonal cycle— as represented 
by the standard deviation (Figure 2b) — spatially covaries with the annual mean 
as well as with climatological expectation. The greatest seasonal variance 
occurs in subtropical regions that experience a marked seasonality in. precipi- 
tation. Central America, Brazil and Africa, for instance, exhibit seasonal 
standard deviations greater than 50 mm. This also is true of Monsoon Asia. 

The suostantial seasonality associated with Mediterranean climates (summer dry) 
also is manifested in lar^e standard deviations of soil moisture along the 

northern coasts of the Mediterranean Sea and the western coast of North America. 

% 

Local and mesos«.nle anomalies are common-^indirectly owing to topographic and 

J # 

vegetatros? influencesyV nonet heless , son moisture variance tends to be inversely 
correlated with the sugnitude of latitude. 

When the soil moisture variance is harmonically decomposed and the har- 
monics are interpreted, the seasonal soil moisture cycle can be understood 
more clearly. Two harmonics were sufficient to explain the majority of the 
variance contained in the seasom.1 soil moisture cycle (Figure 2c). The 
variation left unexplained by the two harmonics (error) is generally between 5 
and 15 on and only at a very few locations does it exceed 20 mm. 

Most of the seasonal variation is explained by the first harmonic whose 
amplitudes frequently approach 100 imn although they (the amplitudes) are highly 
variable in space (Figure 2d). Phase angles associated with the first harmonic 
also are quite variable; once ag:in, owing to the large influence that the 
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local environment exerts on soil moisture. In the mid- to high latitudes of 
the Northern hemisphere, the clockwise rotation of phase angles with latitude 
is partially due to the delay in the snowmelt process which results from reduced 
levels of available energy. That is, since the maximum levels of soil moisture 
often are associated with the infiltration of melt water, the seasonal soil 
moisture peak occurs shortly after the onset of melt. Over much of centra 7 
North America and Asia, this occurs in March or April whereas it occurs in May 
or June in the north of Canada and Siberia. Inland and mountainous environs 
enhance this clockwise rotation, relative to other areas at the same latitude, 
and the effect is evident in the northern Rocky Mountains and Great Plains. 

The contribution of snowmelt to soil moisture in the higher latitudes is 
enhanced by the low magnitudes of potential evapotranspiration early in the 
growing season. 

Equatorward of the seasonal snow line, snowmelt does not influence the 

timing of the soil moisture maxima. Rather, maxima commonly occur at the end 

of extended periods where precipitation exceeds evapotranspiration. Areas 

along the northern shore of the Mediterranean Sea, for instance, experience 

soil moisture maxima in late February and early March, i.e., at the end of the 

rainy season (winter), when potential evapotranspiration is still quite small. 

Deserts such as the Sahara and Gobi exhibit greatly reduced the amplitudes 

associated with the seasonal cycle as they are dry virtually the entire year. 

* . 5 

Subtropical climes, e.g., sub -Sahara^ Africa, Venezuela, and parts of ^outh and 
Southeast Asia, have early-autumn maxima that occur in September or October. 

In many subtropical areas, the seasonal soil moisture cycle has a double 
maximum although the September-October peak is usually dominant. Early-autumn 
maxima result from the northward shift in the ITC2 and accentuated thermal 
convection which responds to the summer accumulation of net irradiance. 
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South of the equator, in subtropical South America, Africa, Borneo and so 
on, soil moisture maxima again coincide with precipitation maxima (Figure 2d); 
that is, they generally occur in February or March. These maxima are nearly 
six months out of phase with their Northern Hemisphere counterparts, and they 
also result from a shift in the ITCZ and thermally induced convection. Further 
south, the mid-latitude deserts, e.g., the Atacama, Kalahari and central 
Australian deserts, are perennially dry and, therefore, have no real seasonality. 
Maxima in southern Chile and Argentina as well as in ^outh Africa and Australia 
occur in August and September^ i.e. , toward the end of winter^ for the same 
reasons that they appear in February and March within similar locales in the 
Northern Hemisphere. In Antarctica, since the soil -arely thaws, no seasonal 
soil moisture cycle is apparent. 

Within most of the seasonally snow covered reaches of the Northern 
Hemisphere, the vector field representing the second soil moisture harmonic is 
quite similar in appearance to the corresponding first harmonic field (Figure 
2c). This suggests that snowmelt -induced maxima are especially well-defined. 
Over eastern North America and northwestern Asia, however, the phase shifts of 
the second harmonic are rotated clockwise with respect to the phase angles of 
the first harmonic in order to account for extended moist periods. In many 
instances, this pattern represents a secondary soil moisture maximum caused by 
summer precipitation which exceeds the environment's ability to evapotranspire 
the moisture. 

Subtropical regions, e.g., in Mexico and sub-Saharan Africa, exhibit 
second harmonic phase shifts which are nearly six months out of phase with the 
first harmonic phase angles. This also means that the second peak (m'^/2 + 6) 

i. 

associated with the second harmonic is in phase with tne first harmonic. Such 
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patterns are a response to the marked precipitation seasonality of these 
regions, i.e., where most of the year is very dry and late summer is quite 


wet — near field capacity/(Win mott et al., 1981 j ). South of the equator, the 
subtropical regions of Africa ec well as much of central South America exhibit 
maxima in May or June which dampen the sub-zero lows in soil moisture predicted 
by the first harmonic in winter (e.g., June, July and August) and lengthen the 
wet-season which extends throughout the remainder of the year. Further south, 
the Southern Hemisphere deserts, e.g., the Atacama, Kalahari and Great Desert 

<y 

of Australia, have a virtually double-'seasonal cycle; that is, of frequency 
two. Eastern Australia and portions of southern central Africa, however, have 
maxima aoout February or March for similar reason s ^ 

In Mozambique, for instance, the f irst^raaximum occurs at the same time of ^ 
year as the phase angle associated with the first harmonic in order to describe 

' r 4 

an accentuated seasonal maximum and dampen the seasonal minimum represented by 
the first harmonic. In eastern Australia, the second harmonic serves essentially 
the same purpose, except that the first harmonic is six months out of phase 
with its Mozambique counterpart. 

Evapo transpiration 

Except in regions that are seasonally water limited for one or more of a 
variety o£ reasons, mean evapotranspiration follows potential evapotranspiration 
which covaries with mean air temperature. For this reason, there is a general 
decrease in Em with the magnitude of .atitude (Figure 3a). Superimposed on 
this large-scale pattern are a number of mesoscale features. 

In the Pacific southwest, southwest South America, central Asia, north 
Africa, southwest Africa and west-central Australia, low precipitation rates 
are responsible for relatively dry soils. This, in turn, provides that nearly 

t 


-18- 


all the available energy is disposed of sensibly and that evapotranspiration 
is less Chan 25 mm/month over most of these zones. Similar average 
evapotranspiration rates occur in the high latitudes but for different reasons; 
that is, cold- climate evapotranspiration is limited by low E° B levels. In 
stark contrast, the wet and seasonally-wet tropics and subtropics exhibit very 
high rates of evapotranspiration. Areas in Brazil, sub-Saharan Africa and 
southeast Asia, for example, evapotranspire well in excess of 100 mm/month. 
Thermally induced precipitation is seasonally or annually plentiful in these 
regions and available energy is relatively high over the entire year. 

Regionally characteristic, annual evapotranspiration rates also are apparent, 
in the seasonally humid east and southern United States. On the whole, 
because of their greater dependence on available energy (E 0 ^), actual annual 
evapotranspiration fields tend to be smoother than corresponding soil moisture 
fields. 

Standard deviations associated with the seasonal evapotranspiration cycles 
rarely exceed 50 mm/month although they are quite variable from region to 
region (Figure 3b). The more pronounced the seasonal temperature and/or 
precipitation cycle the greater the seasonal variance in evapotranspiration. 

When the air temperature falls to near or below freezing during the winter and 
then rises to over 20°C in the summer months and summer precipitation is plentiful, 
a large seasonal variance in evapotranspiration results. Mid-latitude regions 
such as the mid-western and eastern United States, and eastern Asia consequently 
have evapotranspiration standard deviations well over 50 mm/month. Monsoon 
Asia as well as sub-Saharan Africa similarly exhibit standard deviations 
i greater than 50 mm/month. In the first case, this results from coincidence of 
very high levels of E° m and dynamically induced summer precipitation. The 
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concurrence of high siaaer E° m and precipitation also occurs in sub-saharan 
Africa; however, this precipitation is caused thermally rather than 
dynamically. Zones of low variance, on the other hand, occur in areas that 
have little seasonal variance in either precipitation or available energy. The 
rainforest of eastern Brazil, for instance, is warm and moist throughout the 
year whereas the Sahara and Gobi deserts as well as Antarctica are dry over the 
entire year. Tiese regions vary by less than 10 nm/month. Moderately variable 
evapotranspiration cycles characterize much of the rest of the terrestrial 
world . 


Most of the seasonal variability can be described by the first two 
harmonics associated with the evapotranspiration cycle as suggested by the map 
of unexplained variation (Figure 3c). Patterns in the seasonal evapotranspiration 
cycle as expressed by the first harmoni. are coherent and can be explained by 

tutefted ? ; fc 

w a ll - twwmr climatyf principles (Figure 3d). In the mid- and high-latitudes, 

most of the seasonal variance in Em is explained by the first harmonic since 

E° m follows the sun and, therefore, is inversely proportional to the magnitude 

of the latitude. It follows that the mid- and high-latitude evapotranspiration 

maxima should occur during the high-sun season— which they do. On the more 

humid, eastern portions of the continents, a lag (approximately one month) in 

the evipotranspiration maxima— behind the solar maximum— corresponds to the 

familar peak in air temperature which also occurs after the solar maximum. In 

lower latitudes, for example central America, evapotranspiration 

maxima occur even later in the suumer ^e.g., in August and September) since the 

amplitude of the seasonal solar cycle is lessened and the accumulation of net 

irradiance continues through- out the summer This same effect also is apparent 

S 

in sub-Saharan Africa and jtouth Asia as well as in the Southern Hemisphere. 


Note the clockwise rotation of phase angles in Brazil, southern Africa and 
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Australia as the equator is approached. Apparent anomalies arise in 9teppe ov 
Mediterranean areas (e.g., in northern Libya) since maxima largely depend on 
moisture, rather than temperature, which may precipitate in but a very few of 
the winter months. Spring maxima east of the Caspian Sea also result ‘from a 
shortage of precipitation and subsequently of soil moisture. Evapo transpiration 
maxima that do not eovary with the seasonal maxima of available energy (e.g., 
air temperature) usually have large soil moisture deficits during part of the 
year or a nearly constant available energy supply. 

Especially in the higher latitudes of the Northern Hemisphere, interpreta- 
tion of the second evapotranspiration harmonic siay seem complex owing to 
frequent reversals in the phase shift (Figure 3e). It should be remembered, 
however, that harmonics of frequency two have two maxima separated by six 
months and, therefore, neighboring grid locations that seem to be out of phase 
by five to six months are really only out of phase by less than one month. 

Over much of uorthern Europe, Asia and North America, therefore, secondary 
evapotranspiration maxima are evident in June and July. These summer maxima 
describe an increase in the intensity of the summer evapotranspiration season 
with latitude — relative to negligible winter evapotranspiration rates. Corres- 
ponding maximum phase shifts that occur in December and January (i.e., six 
months out of phase with the summer maxima) when added to the first harmonic 
serve to dampen the erroneous, negative winter evapotranspiration rates 
predicted by the first harmonic. Similar patterns also are apparent in 
Antarctica. As one moves south from the Northern Hemisphere high latitudes 
into the mid-latitudes, a counterclockwise rotation of the phase shift is 
apparent, particularly in the western United States, Europe and western Asia 
(Figure 3e). This shift represents a lengthening of the evapotranspiration 
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season owing Co Che increases in available energy wich Che decrease in lacicude 
as well as Co Che marine influence on Che vescern reaches of Che concinencs. 

On Che easCern porCions of North America and Asia, Che lengCh of Che warm season 

.S 

is affecCed liCCle by adjacent oceans. Similar paCCerns exist in Che Southern 

H 

hemisphere mid-lati' ides, although they are not so well-developed as in the 
norCh because of Che relatively limited landmass of Che mid-latitude Southern 
Hemisphere . 

In Che subtropics, a substantial increase of variance wich longitude, in 
both Che phase shifts and amplitudes associated wich the second harmonic, is 
apparent (Figure 3e). Over much of north Africa, for example, Che seasonal 
cycle is virtually not present whereas it is exaggerated by Che Southeast Asian 
monsoon. The subtropics of the Southern hemisphere by contrast exhibit a 
relatively coherent extension of the evapotranspiration season well into spring 
and fall. Evapotranspiraiion in these areas, of course, occurs throughout most 
of the year. 

In much of the tropics, the amplitudes of the second harmonics are smaller 
than those of their subtropical and mid-latitude counterparts which represent^ 
a dampening of the seasonal cycle (Figure 3e). The phase «hifts become more 
variable in an effort to describe a perennial evapotranspiration regime whose 
maxima often are not well-defined and may be subject more to temporally variable 
local and mesoscale influences than to the large-scale, solar-fotced seasonal < 

A 

cycle. Equatorial South America illustrates this variability, while equatorial 
Africa portrays a greater degree of coherence because of its lesser topographic 
variability. 
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CONCLUDING REMARKS 

Using recently available monthly air temperature and precipitation data 
for the vorld and an empirically-baaed water budget algorithm, the large-scale 
spatial and seasonal variability associated with terrestrial snow cover, soil 
moisture and evapotranspiration have been estimated and inductively evaluated. 

Even though ve regard these findings as a zeroth-order approximation — one of a 
very few available, these data provide both empirically-abased climate fields <^— 
to which general circulation results can be compared and realistic initial 
fields that may be used in GCM seasonal cycle experiments. 

Of the three climate fields presented, the seasonal snow cover cycle 
exhibits the greatest statistical regularity, particularly in the higher latitudes 
of the Northern Hemisphere. In the mid-latitudes, on the other hand, the 
seasonal snow cover as well as its albedo and water equivalent can be highly 
variable, both spatially and temporally. As snow has a high albedo and the 
water entrained in the snowpack ie a potentially significant source of soil 
moisture (especially during the melt season which closely precedes the growing 
season) these results suggest that the interrelationships between the seasonal 
snow cycle and large-scale atmosphere circulation need to be well-represented 
within GCMs . For a number of reasons, it is not sufficient to designate snow 
cover as either present or absent. The seasonally transient water contained in 
the snowpack, fer example, has an Important influence on the surface energy 
balance through latent flux which derives from melt water; that is, snow cover 
has climatic significance beyond its obvious effect on regional albedo. The 
magnitude of the sno:-* water store, in fact, can be as large as or larger than 
the soil moisture store. Since the influence of snow and snow water equivalent 
on climate is affected by topographic and vegetational factors which have not 
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been considered here, further research into the large- scale role which the 

terrestrial snow cycle plays in the seasonal climate cycle is required. It is 

known, for example, that the regional albedo of a snow covered forest is 

substantially less than a snow covered grassland (Leonard and Eschner, 1968) 

which, of course, directly affects the surface energy balance. The differential 

melt and subsequent soil moisture responses of seasonally snow covered forests 

and grasslands as well as ensuing effects on the l&rge scale circulation of the 

atmosphere are not so well-known. It is likely that more detailed representations 

of the seasonal snow cycle will unveil more realistic but even more variable 

snow cover and snow water equivalent fields. 

Terrestrial soil moisture is much more temporally and spatially variable 

than snow cover or snow water equivalent because of its greater dependence upon 

the ensemble of local environmental factors. A number of large-scale patterns 

+w ** v 

are evident, however, can be plained by accepted climate principles. The 

r 

mid- to high latitudes of the Northern Hemisphere, for example, reach peak soil 
moistures after the snow accumulation season during the spring thaw while soil 
moisture in near monothermal tropics follows the seasonal rainfall cycle. Even 
though our maps suggest that the seasonal soil moisture cycle is highly variable, 
it is probable that more realistically-derived fields that recount for soil, 
vegetational and physiographic variables will exhibit even greater variability. 
Future representations of terrestrial soil moisture, as a consequence, should 
include more of the salient near-surface environmental and hydrologic variables. 
Seasonal evapotranspiration has less spacio-temporal variability than soil 


moisture but more than snow cover. The former is principally due to the observation 
that evapotranspiration is relatively insensitive to changes in soil moisture 
until soil moisture level falls well-below one-half the field capacity (e.g., 
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see Mather, 1974; 106: Nappo, 1975). Once again, our characterization may be 
oversimplied, since, the vegetation structure and morphology as veil as 
physiographic and pedological considerations are absent. 

Because the spacio-temporal variability of evapotranspiration has key 
importance for the general circulation of the atmosphere (Mintz, 1982), the 
relationships between the terrestrial biosphere and the near-surface atmosphere 
should be modeled in greater detail (Mintz et al., 1983). At the same time, 
empirically-derived evapotranspiration climatologies must be improved, most 
likely through remote sensing means. Nevertheless, the physical expectation of 
the macro- and mesoscale features of the terrestrial seasonal evapotranspiration 
cycle (e.g., evapotranspiration follows available energy except in areas that 
are seasonally or annually water-limited) seems to be reasonably well- 
represented in this climatology of evapotranspiration. 
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FIGURE CAPTIONS 


Figure l Baaed upon seasonally representative, mid-monthly estimates of snow 
water equivalent at the nodes of a 4° of latitude by 5° of longitude 
lattice, map 8 (a) and (b) are isoline representations of the means 
and standard deviations, respectively, associated with the seasonal 
cycles at the grid points. The variances contained within tho grid 
point seasonal cycles were harmonically decomposed and, for all grid 
points, the roots of the variances left unexplained by the first two 
harmonics (i.e., the errors) were contoured and are presented in map 
(c). Haps (d) and (e) give the spatial distributions of the amplitudes 
(lengths of the arrowed vectors) and phase shifts (direction of the 
arrow heads) associated with the first and second harmonics, respectively. 


Figure 2 Based upon seasonally representative, mid-monthly estimates of 

_f 0 1 1 . Q 

♦vapor rang pirate*, o tt at the nodes of a 4° of latitude b}' 5° of 

A. 

longitude lattice, maps (a) and (b) are isoline representations of 
the means and standard deviations, respectively^ associated with the < 
seasonal cycles at the grid points. The variances contained within 
the grid point seasonal cycles were harmonically decomposed and, 
for all grid points, the roots of the variances left unexplained by 
the first two harmonics (i.e., the errors) were contoured and are 
oresented in map (c). Maps (d) and (e) give the spatial distributions 
of the amplitudes (lengths of the arrowed vectors) and phase shifts 
(direction of the arrow heads) associated with the first and second 
harmonics, respectively. 
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Figure 3 Based upon seasonally representative, mid-monthly estimates of 

-ee rt me i s cu rerat the nodes of a 4° of latitude by 5° of longitude <£T~ 
lattice, maps (a) and (b) are isoline representations of the means 
and standard deviations, respectively, associated with the seasonal 
cycles at the grid points. The variances contained within the grid 
point seasonal cycles were harmonically decomposed and, for all grid 
points, the roots of the variances left unexplained by the first two 
harmonics (i.e., the errors) were contoured and are presented in 
map (c). Maps (d) and (e) give the spatial distributions of the 
amplitudes (lengths of the arrowed vectors) and phase shifts 
(direction of the arrow heads) associated with the first and second 
harmonics, respectively. 


FIGURE TITLES 


Figure 1 


Figure 2 


Figure 3 


Seasonal Water Equivalent of the Snowpack 

a) Annual Average (mm) 

b) Seasonal Standard Deviation (mm) 

c) Root of the Variance Unexplained 
by the First Tvo Harmonics (mm) 

d) Amplitude (mm) and Phase Angle 
of the First Harmonic 

e) Amplitude (mm) and Phase Shift 
of the Second Harmonic 


Seasonal Soil Moisture 

a) Annual Average (mm) 

b) Seasonal Standard Deviation (mm) 

c) Root of the Variance Unexplained 
by the First Two Harmonics (mm) 

d) Amplitude (mm) and Phase Angle 
of the First Harmonic 

e) Amplitude (mm) and Phase Shift 
of the Second Harmonic 


Seasonal Evapotranspiration 

a) Annual Average (mm/month) 

b) Seasonal Standard Deviation (mm/month) 

c) Root of the Variance Unexplained by the First 
Two Harmonics (mn/month) 

a) Amplitude (mm/month) and Phase Angle of the 
First Harmonic 

e) Amplitude (mm/conth) and Phase Shift of the 
Second Harmonic 
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APPENDIX 3 


Evaluation of Thomthwaite 's Model Through 
Comparisons with Lettau's Model. 



EVALUATION OF THORNTHWAITE ' S MODEL THROUGH COMPARISONS WITH LETTAl' ' S MODEL 


Phase three of the project called for an evaluation or the Thornthwaite 
model of actual eve po transpiration under climatic conditions where it is felt 
to provide questionable results. In addition, the model is often validated 
by comparing measured stream runoff with that predicted by the model. 
Discrepancies between the time series of calculated and predicted runoff are 
reduced by adjusting the soilwater holding capacity or the 3 function- Although 
one might assume that the resulting soil moisture time series is correct, 
this procedure may produce compensating errors in evapotranspiration and soil 
moisture storage. It is therefore desirable to validate Thornthwaite ' s model 
determinations of actual evapotranspiration and soil moisture in a variety of 
climates. 

For example, for the Maurice River basin in southern New Jersey, during 
the years 1950-1956, the total precipitation was 8204 mm (1172 mm/yr) and the 
total measured runoff was 3398 mm (per unit area). Assuming that the net 
change in water storage over these seven years is negligable and that there 
is no subterranian discharge from the basin, the total evapotranspiration for 
the period is 4806 mm. On the other hand, the total evapotranspiracion cal- 
culated by the Thornthwaite model for the period depends on the assumed value 
cf the model soil moisture storage. With an assumed value of storage of 25 
mm, the accumulated evaporation is 4337 ram. Increasing the assumed storage 
to 225 mm increases the accumulated evapotranspiration to 4781 ram. The smaller 
value is 10% below the difference between precipitation and runoff- 

In the case of a dryer climate, for the experimental agricultural water- 
shed #111 in Chikasha, Oklahoma during the years 1962-1972 there was an 
accumulated precipitation of 6 1 5 1 mm (599 tran/yr), and a runoff of 2f>9 mm 
giving an accumulated evapotranspiration of 5951 mn. . Accumulated evapotrans- 



piration for this site according to the Thornthwaite method varies from 5662 
mm when a storage of 125 is assumed to 5978 mm for a storage of 250 mm. The 
former is 51 belov the difference of the accumulated precipitation and runoff. 

In addition, it is known that the Thornthwaite method underestimates 
cold season evapotranspiration. Therefore, if the accumulated amount is ad- 
justed to agree with accumulated precipitation less accumulated runoff, it 
must over estimate warm season evapotranspiration by an equivalent amount;. 

The most frequent means for checking model performance has been to com- 
pare the model with lysimeter observations of actual evapotranspiration. 

Such an approach is not satisfactory for a study involving large areas of 
landscape because of the impossibility of relating lysimeter measurements, 
usually of a single agricultural crop, to the evaporation from deep rooted, 
often forested, vegetation. Uncertainty is increased by the heterogeneous 
character of plant soil and terrain features found in extended real landscapes. 

A possible means for checking the Thornthwaite model has been suggested 
by an approach developed by Lettau (1969) and Lettau (personal communication) 
and developed by Lettau and Baradas (1973), Lettau, Lettau and Molion (1979). 

The aim is to explore the possibility that this model can be calibrated 
against actual watershed measurements to derive the seasonal time series of 
actual evapotranspiration and soil moisture. Lettau' s Model is developed 
from a water budget equation 

6m 

P = N + E + (5) 

where P, N, and E are the fluxes of precipitation, runoff and evapotranspiration 
per unit area per unit of time ( t) (day, month, year, etc). Moisture, m, is 
the exchangeable moisture stored anywhere in the soil, whether in the root 
zone or not, so long as it is subject to eventual removal from the stream basin. 

Equation (5) is parameterized by dividing runoff and evapotranspiration 
into components which occur in the same t 1 me interval in which the rainfall 


3 


occurs (called immediate runoff and evapotranspiration and designated by 
single primes) and runoff and evapotranspiration which is derived from preci- 
pitation falling in previous time intervals (called delayed runoff and evapo- 
transpiration and signified by double primes). Thus 

N ■ N’ + N" 

E » E’ + E" 


The parameterized form of equation (5) is 

? 


6m 


P - N' - E' «• P* + N" + E" + (6) 

The respective terms are defined in terms of the six hydrologic parameters 
N’ 


n p (P - Pn); N" = vN™ 


E' » e p (P - Pe); E" = ve 


,m 


(7) 


where n p and e p are dimensionless coefficients between 0 and 1 which dete mine 
the fraction of the precipitation in excess of the thresholds Pn and Pe that 
actually run off, and vN and vE are coefficients with units of reciprocal time 
which describe how frequently the stored soil water will be exchanged during 
the time interval At. The value of vN and vE will be greater than 0 and usually 
but not necessarily less than 1. 

Lettau derives a differential equation for soil moisture from these ideas 
which, when integrated in a stepwise forward manner, yields 

m. , - (P - N' - E ' ) / v. + (m. - (P - N 1 - E;)/v.) exp (-v.At) (8) 

l+l j i J J 

where nu is the stored exchangeable soil water of the beginning of the ^th time 
interval, 


v. “ vN. + vE . 

J J J 

From this time series of exchangeable soil moisture one may then calculate 
evapotranspiration and runoff according to equation 7. 

The parameters may be modified in various ways to improve the representation 
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of the physical processes. vN may be a function of soil moisture. The two 

evaporation coefficients e^ and vE are likely to be functions of absorbed 

solar radiation. This introduces the need to include albedo, which may 

change with seasons and/or be a function of soil moisture. Further n or Pn 

P 

may be made functions of soil moisture. 

The Lettau model may be used as an independent test of the Thomthwaite 
water budget calculations of evapotranspiration and soil moisture. However, 
this will require knowing the six parameters, as functions of time if necessary, 
that describe the hydrologic behavior of stream basins for which good measure- 
ments of runoff, precipitation, and absorbed solar-radiation are known. It 
also requires developing an algorithm for determining the set of six para- 
meters that best fit the hydrologic regime of a particular stream basin. 

Finding a coherent set of parameter values which will allow the model to 
describe the hydrolic behavior of a particular river basin is called cali- 
brating the model. This could be accomplished unambiguously if we knew the 
four time series of immediate and delayed runoff and immediate and delayed 
evapotranspiration. The time series of soil moisture storage could be sub- 
stituted for one of the delayed process series. If the time series of precipi- 
tation, absorbed solar radiation, and runoff are known or estimable for a 
particular basin, it is possible to obtain sets of parameters whic. permit 
calculating the unknown time series. (Lettau and Baradas , 1973; Grangier, 

1976, 1983). 

The time series of immediate and delayed runoff (N 1 and N") can be obtained 
from a stream gage record of discharge. Several approximation procedures have 
been used. However, for this work a more accurate procedure has been developed. 
The new computerized procedure uses the daily discharge and precipitation record 



5 


over a period of years to calculate the immediate runoff following a single 
or group of precipitation events which is then aggregated into monthly sums. 

The new method produces a reliable partitioning for the data examined in this 
study. Having obtained N' and N" by differencing, the parameters describing 
them can later be obtained when the rest of the analysis is completed. 

The two evaporative tine series remain unknown, as do the series of soil 
moisture ,m£ and A.m£ . However, if we assume 0 we can obtain the total- 

evapotranspiration from the difference of the precipitation and runoff sums. 
This amount must be partitioned between immediate and delayed processes. 

Lettau and Baradas suggest a method of estimating a physically plausable guess 
for the immediate evaporicity, e* , which is used here. They then obtain time 
series of E" and m£ by a successive approximation routine which uses the time 
series of N" to improve the estimates. The routine converges very quickly. 

In '.he case of their tropical data, the resulting values of E" and m£ are" 
plausable. The retentivity, v, is found to be a function of m£ . When this 
method is applied to the data from New Jersey, the result is not satisfactory 
because the calculated soil moisture and m£ series show no seasonality. It 
appears that their procedure works when the seasonality is in the precipitation 
regime but noc when the seasonality is in the solar radiation (temperature) 
regime, as in a temperate climate. 

An alternative separation procedure is required which will produce the 
observed evapora.tion and soil moisture seasonal patterns. Such a procedure 
cannot be based solely on the seasonal march of absorbed solar radiation be- 
cause that change is not large enough (here about 3.6:1) to account for the 
seasonal variation in evapotranspiration (5:1). It appears that the phono- 
logical Cj -:le of the vegetation actively regulates the evapotranspiration. 
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Several procedures to obtain E" were tried, none really satisfactory. 

The one presented here is preferred to the others because of its similarity 
to Lettau's approach which is based on the premise that a solution has been 
found if the chosen parameter values reveal a regular functional relationship 
between model soil moisture and the delayed process parameters. 

The approach of Lettau and Baradas was modified to assert that the delayed 
evapc transpiration be proportional to the normalized solar radiation, the 
soil moisture, and a new factor ct which describes the relative effectiveness 
of the vegetation in transpiring. A value of 1 is applied for the months of 
full vegetation development, June, July, August and September. A value of .1 
is applied in December, January and February; .2 for March and November; and 
.7 for May and October. Alpha (a) is applied only to the delayed evaporation 
during the calibration procedure. It implies the assumption that vegetation 
is responsible for the delayed evapotranspiration, while che immediate eva- 
potranspiration proceeds independently of the stage of the vegetation. 
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Note that the soil moisture series may be calculated numerically , 6y forward 
integration, according to Lettau's solution (eq. 8) 
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The iteration is continued, with the observed delayed runoff being used to 
revise the time series of the partitioning parameter u* and until the 
calculated values of N M equal the observed values of N". Closure occurs 
within three or four rounds. The results of using this procedure on the 
data for the Maurice River in New Jersey will illustrate the type of result 
which can be achieved. 

The Maurice River in southern New Jersey has a drainage area of about 
100 square miles, of which 48% is agricultural land, 29% oak-pine forest, 

20% hardwood swamp forest, and 1% each broad leaf upland forest, urban land, 
and open water. The soils for 85% of the area are graveley and sandy, and 
are well drained, lne terrain is gently sloping coastal plane. A weighted 
average for the water holding capacity of the soils lies between a low of .1 
and a high of .17 per cm of top soil. The river is dammed at Norma, where 
the stream gage is 1 ted. A lake and marshes extend upstream from the dam. 

Daily precipitation for Clayton, located in the center of the basin, was 
used. (Data from other surrounding stations were used in monthly analysis 
only). The daily gaged runoff for the years 1949-1962 was separated into 
immediate and delayed runoff by the routine described previously. Monthly 
solar radiation was obtained from the Seabrook station, about 12 miles distant, 
for most of the period. Missing months were estimated from a regression on 
an average of New York City and Washington, the coefficient of determination 
for which against Seabrook is .982. Albedo for the basin was obtained as a 
weighted average of values obtained from the literature for such land uses 
and seasons. The principle reference used was Rung, Bryson, and Lenschow 
(.1964). The region has relatively mild winters with average monthly tempera- 
tures only infrequently belcw freezing. Snow rarely lingers longer than five 
days. The soil may freeze if snow does not accompany the cold. 
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The years 1949-1962 were divided into years of below average pre ipi- 
tation (dry) and years of above average precipitation (wet) for separate 
processing. The data and analysis for wet, dry and all years are given 
below. (Tables 1 and 2) Figures 1 and 2 are the plots of the hydrologic 
data used in attempting to calibrate the model. There appears no reason 
for uncertainty in the separation of total runoff into immediate and delayed 
runoff. ’Figure 3 and 4 are the graphic output of the interactive microcomputer 
program which is my implementation of a calibrating algorithm. For the case 
illustrated, the total flushing frequency v was initially set to .25 for 
both wet and dry years and e* (HE) is .15. Figure 3 is a plot of the monthly 
values of Vjj * N"obs/' J i obtained after convergence (W for wet years and D 
for dry years). Notably, v N appears to be independent of soil moisture. 

This result is a consequence of the procedural assumption for separating out E". 
Further, there is considerable scatter of point3 which probably is also an 
artifact of forcing all of the soil moisture dependence onto E". Figure 4 is 

F i 

a plot of normalized E" against the soil moisture where normalized E" = E"/ct — . 
In this case the linear plot is forced by the procedure and does not reveal 
anything about the quality of the fit produced by the selected parameters. 
Ideally, these tow plots would reveal the best fitting set of parameter values 
by showing a tendency toward smooth functions as the proper parameter values 
are chosen. In these runs there were no clear trends toward functional re- 
lations, so it is not possible to choose proper values for e* and v extent to 
assert that they each probably lie within the range of .15 to .3. 

Figure 5 shows a plot of the calculated value of the total evapohranspi- 
ration in two cases, that where v * .15 and e* = .15 d where v - .25 and 
e* * .25. Values of actual evapotranspira., ion cal ated according to the 

Thornthwaite model for the years 1950-1956, with SLorage 225 mm are also 
given for comparison. 
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An examination of the three evapotranspiration plots in figure 5 reveals 
first that the two Lettau models show very little difference in the amount 
and timing of the total evapotranspiration. The model with the larger 
flushing rate and larger immediate evaporicity shows a small tendency for 
larger evapotranspiration early and late in the year and a slightly lower 
peak in June. Both show the strong dependence, on absorbed solar radiation, 
which reaches its maximum value in June. Total evapotranspiration then falls 
rather rapidly as the soil moisture is drawn down and the absorbed solar radia- 
tion diminishes. The Thomthwaite model actual evapotranspiration shows con- 
siderably less water loss in the winter and spring, peaks in July and continues 
to show higher evapotranspiration than the Lettau model until December. The 
timing of the water loss is shifted according to this models dependence on 
air temperature. 

There is no clear basis for choosing between these two models. The factor 
a used in the Lettau model could be adjusted to change the seasonal timing 
toward that of the Thomthwaite model. However, the values actually used were 
chosen according to estimates of the state of vegetation development in each 
month (see Tables 1 and 2), and should be changed only on the basis of better 
objective information about the relative seasonal transpirational effective- 
ness of the vegetation. 

There is a need to devise a better procedure for obtaining the Lettau 
delayed evapotranspiration in climates such as this where the seasonality is 
driven by the available energy and not by the precipitation. Until then a 
precise basin calibration for seasonal cycles cannot be done. 
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